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Abstract 

The paper deals with Delaunay Triangulations (DT) in E d space. This classic 
computational geometry problem is studied from the point of view of the efficiency, 
extendibility to any dimensionality, and ease of implementation. 

A new solution to DT is proposed, based on an original interpretation of the well- 
known Divide and Conquer paradigm. One of the main characteristics of this new 
algorithm is its generality: it can be simply extended to triangulate point sets in 
any dimension. The technique adopted is very efficient and presents a subquadratic 
behaviour in real applications in E 3 , although its computational complexity does 
not improve the theoretical bounds reported in the literature. An evaluation of the 
performance on a number of datasets is reported, together with a comparison with 
other DT algorithms. 
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1 Introduction 


Triangulation is one of the main topics in computational geometry and it is commonly used in a 
large set of applications, such as computer graphics, scientific visualization, robotics, computer 
vision and image synthesis, as well as in mathematical and natural science. Given a point set P, 
the Delaunay Triangulation (DT) is a particular triangulation, built on the points in P , which 
satisfies the empty circum-circle property: the circum-circle (-sphere in E 3 or -hypersphere in 
E d ) of each simplicial cell in the triangulation does not contain any input point p £ P. Many 
algorithms have been proposed for the DT of a set of sites in E 2 , E 3 or E d , and most of them 
are reviewed in [1, 2]. 

Unfortunately there has been little research into implementations and performance eval¬ 
uations of Delaunay triangulators. Few papers report evaluations of real implementations or 
give experimental comparisons of different algorithms. Worst case time complexities are gen¬ 
erally given, but such analyses, from the point of view of the application programmer, are 
not always sufficient to make the correct decisions. In fact, theoretically better algorithms can 
sometimes be outperformed by more naive methods; the theoretical asymptotic worst case com¬ 
plexity sometimes fails to consider the optimization techniques that can be applied to reduce 
the expected complexity. 

A new divide & conquer DT algorithm is proposed in this paper. The algorithm gives a 
general and simple solution to DT in E d space and makes use of accelerating techniques which 
are specific to computer graphics. 

The paper is organized as follows. Definitions and a taxonomy of Delaunay triangulation 
algorithms are presented in Section 2. The proposed algorithm is described in detail in Section 3, 
together with some optimization techniques. The performances of the proposed solution are 
evaluated on a number of datasets and compared with other solutions in Section 4. Conclusions 
are drawn in the last section. 


2 Delaunay Triangulation 

Given a point set P in E d , a k-simplex, with k < d, is defined as the convex combination of k + 1 
affinely independent points in P, called vertices of the simplex (e.g., a triangle is a 2-simplex 
and a tetrahedron is a 3-simplex). An s-face of a simplex is the convex combination of a subset 
of s + 1 vertices of the simplex (i.e., a 2-face is a triangular facet, 1-face is an edge, 0-face is a 
vertex). 
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A triangulation £ defined on a point set P in E d space is the set of d-simplices such that: 

1. a point p in E d is a vertex of a simplex a in £ iff p £ P; 

2. the intersection of two simplices in £ is either empty or a common face. 

3. the set £ is maximal: there does not exist any simplex cr that can be added to £ without 
violating the previous rules; 

A triangulation £ is a Delaunay Triangulation iff the hypersphere circumscribing each 
simplex does not contain any point of the set P [3, 4]. The Delaunay triangulation of a given 
point set P is unique if there do not exist in P d + 2 points lying on the same hypersphere. 
Such cases, also known as degeneracies, can be managed by using local perturbation schemes 
[5], 

The duality between DTs and Voronoi diagrams is well known [4] and therefore algorithms 
are given for the construction of DT from Voronoi diagrams. However, direct construction 
methods are generally more efficient because the Voronoi diagram does not need to be computed 
and stored. Direct DT algorithms [1] can be classified as follows: 

• local improvement - starting with an arbitrary triangulation, these algorithms locally 
modify the faces of pairs of adjacent simplices according to the circum-sphere criterion; 

• on-line (or incremental insertion) - starting with a simplex which contains the convex 
hull of the point set, these algorithms insert the points in P one at a time: the simplex 
containing the currently added point is partitioned by inserting it as a new vertex. The 
circum-sphere criterion is tested on all the simplices adjacent to the new ones, recursively, 
and, if necessary, their faces are flipped; 

• incremental construction - the DT is constructed by successively building simplices whose 
circum-hyperspheres contain no points in P; 

• higher dimensional embedding - these algorithms transform the points into the E d+1 space 
and then compute the convex hull of the transformed points; the DT is obtained by simply 
projecting the convex hull into E d ; for a comparison of the different approaches see [6]; 

• divide & conguer (D&C) - this is based on the recursive partition and local triangulation 
of the point set, and then on a merging phase where the resulting triangulations are joined. 
Current algorithms are not generalized to E d space, but limited to E 2 space alone. 








On-line methods [7] hold the lower worst case time complexity, 0(n log n + n^-l) [8]. 
Moreover, these methods in their naive implementation are simple to program and can be 
generalized to manage point sets in E d space. 

D&C solutions hold in E 2 the same complexity as on-line methods, but a general D&C 
E d (d > 2) solution has not been proposed yet. The main problem here is the design of the 
merging phase. Due to the explicit ordering of the edges incident in a vertex (Figure 1), the 
merging phase is simple in E 2 [9], but hard to design in E d where this ordering is not given. 

The algorithm proposed in this paper bypasses this problem by reversing the order between 
the solutions of sub-problems and the merging phase. The classical D&C algorithms recursively 
subdivide the input points, construct two partial DTs and then merge them. Our solution is 
based on a more complex division phase, in which the input dataset P is split into Pi and P 2 , 
and a section of the DT is immediately built. This partial triangulation allows the algorithm 
to recursively triangulate the two point sets Pi and P 2 , taking into account the border of the 
partial triangulation and avoiding the need for a further merging phase. A “merging” simplex 
set is thus built before the subproblems are solved: we partition the problem solution, instead of 
its instance. The partial triangulation can be built very simply using a constructive rule similar 
to Me Lain’s in its incremental construction approach [10]. This means we can specify a general 
E d D&C Delaunay triangulator. Its simple structure permits an efficient implementation using 
some well known optimization techniques. 

3 The DeWall Algorithm 

A new algorithm for the DT of a point set P in E d is presented in this section. The algorithm 
is based on the D&C paradigm, but this paradigm is applied in a different way with respect 
previous DT algorithms [9] [11]. The general structure of D&C algorithms is: divide the input 
data into subset Pi and P 2 ; recursively solve on Pi and P 2 ; and merge the partial results Si 
and S2 to build solution S. 

In the case of triangulations, the input point set P can easily be divided using a cutting 
plane such that the two associated halfspaces contain two point sets Pi and P 2 of comparable 
cardinality. The problem is how to implement the merging phase, i.e. how to build the union 
of the two solutions Si and S 2 . This union requires the triangulation of the space separating Si 
and S 2 , and generally also requires a number of local modifications to Si and S 2 . As previously 
stated, this problem has been efficiently solved for the E 2 case [9] [11], but not for the general 
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Our approach to D&C is slightly different. Instead of merging partial results, we apply a 
more complex dividing phase which partitions the point set and builds, as first step, the merging 
triangulation. The algorithm is then recursively applied to triangulate the two subsets of the 
input dataset P. 

The splitting plane a separates the point set P into two subsets Pi and P2. Analogously, 
the splitting plane a divides a triangulation £ into three disjoint subsets: the simplices that 
are intersected by the plane, which we call the simplex wall £„, and the two sets of simplices 
£1 and £2 that are completely contained in the two halfspaces defined by a (Figure 2). £„ can 
be choosen as a valid merging triangulation: (a) each <7 £ £„ is also in £ and (b) subtracting 
£„ from £ generates two disconnected simplicial complexes £1 and £ 2 . 

The DeWall (Delaunay Wall) algorithm, specified in pseudo-pascal in Figure 3, consists 
of the following steps: 

• select the dividing plane a, split P into the two subsets Pi and P 2 and construct £„; 

• starting from £„, recursively apply DeWall on Pi and P 2 to build £1 and £ 2 ; 

• return the union of £„, £1 and £ 2 . 

The technique used to build the simplex wall £„ is a slight variation on an incremental con¬ 
struction algorithm; it is described in the next section. 

3.1 Incremental construction of the simplex wall 

The simplex wall can be simply computed by using an incremental construction approach: a 
starting simplex is individuated and then £ is built by adding a new simplex at each step 
and without having to modify the current triangulation. This technique for DT was originally 
proposed in E 2 by Me Lain [10], and then applied by Dobkin and Laszlo [12] for E 3 subdivisions. 

The incremental construction approach can be easily generalized to E d triangulations: for 
each (d— l)-face /, which does not lie on the Convex Hull(P), there are exactly two simplices a 1 
and < 7 2 in £, such that oq and <7 2 share the (d — l)-face /. The algorithm starts by constructing 
an initial simplex cq; then, it processes all of the (d — l)-faces of cq: the simplex adjacent 
to each of them (if it exists, i.e. the face does not belong to the Convex Hull of P) is built 
and added to the current list of simplices in £. All of the new (d — l)-faces of each new 
simplex are used to update a data structure, here called Active Face List (AFL). Update of the 
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AFL is as follows: if a new face is already contained in AFL, then it is removed from AFL; 
otherwise, it is inserted in AFL because its adjacent simplex has not yet been built. The process 
continues iteratively (extract a face / from AFL, build the simplex cr adjacent to /, update the 
AFL with the [d — l)-faces of cr, and then again extract another face from AFL) until AFL is 
empty. In the implementation of the AFL data structure, the efficiency of the most common 
operations (Insert, Extract, Delete, Member) has to be guaranteed. Our implementation of the 
AFL data structure is based on hash indexing, making it possible to manage AFL in nearly 
constant time (an average of 1.15 - 1.5 accesses to the hash table were measured with the current 
implementation to solve a query). 

Given this general incremental construction algorithm, we only need to specialize it for the 
construction of E„. In particular we have to detail: (a) how to build the initial simplex (the 
MakeFirstSimplex function), (b) how to build the simplex adjacent to a face / (the MakeSim¬ 
plex function), and(c) how this construction process can be limited to the simplices in E a . 

Construction of the first simplex 

The function MakeFirstSimplex produces a Delaunay d-simplex which is intersected by the 
plane a, in order to start from this simplex the incremental construction of the simplex wall 

MakeFirstSimplex selects the point pi £ P nearest to the plane a. It then selects a second 
point P2 such that P2 is the nearest point to p\ on the other side of a. Then, it searches the 
point p 3 such that the circum-circle around the 1-face (pi,p 2 ) and the point p 3 has the minimum 
radius; (pi,P2,P3) is therefore a 2-face of E. The process continues until the required d-simplex 
is built. 

Construction of the generic simplex 

Given a face /, the function MakeSimplex builds the adjacent simplex by applying the DT 
definition. For each point p £ P, MakeSimplex computes the radius of the hypersphere which 
circumscribes p and the face /. We choose the point p which, generally speaking, minimizes 
this radius to build the simplex adjacent to /. 

MakeSimplex selects the point p which minimizes the function dd (Delaunay distance): 

j r if c E Halfspace(f,p) 
dd(f,p) = < 

I — r otherwise 


with r and c the radius and the center of the circumsphere around / and p; given the plane on 






which / lies, Halfspace(f,p) returns the halfspace which contains the new tetrahedra. 

Let us introduce the following example to illustrate the definition above. Let us assume that 
there exists a subset of points Q, which are contained in Halfspace(f,p) and are located on a 
straight line which intersects the face /. If these points are processed in order of decreasing 
distances from /, and the centers of the circumspheres are computed, we can observe that these 
circumsphere radii will decrease until we get the first point rp £ Q whose circumsphere center is 
located in the opposite halfspace (the one which do not contains the points in Q). Then, for all 
successive points qu,k > i, the radii will start to increase. The dd distance defined previously 
takes into account this decreasing-increasing behaviour of the circumsphere radius. 

The analysis of the points p £ P is limited to the points which lie in the outer halfspace with 
respect to face / (i.e. the halfspace which does not contain the previously generated simplex 
that originates face /). 

The outer halfspace associated with / contains no point of P iff face f is part of the Convex 
Hull of P (the faces on the Convex Hull are the only faces that belong to just one simplex in 
the triangulation). In this case the algorithm correctly returns no adjacent simplex and, in this 
case only, MakeSimplex returns null. 

A simple solution to reduce the cost of MakeSimplex function is to take into account, for 
each point p in P, the current triangulation progress status. As soon as all of the simplices 
incident in p have been built, p may be removed from P and it will no longer be tested in the 
further invocations of MakeSimplex. The control on the number of incident simplices was im¬ 
plemented with a counter associated with each vertex p, increased each time a new face incident 
in p is built and decreased for each invocation of MakeSimplex on an incident face; as soon as 
the counter returns zero, p may be deleted from P. 

Construction of simplices in S a alone 

A slight modification to the canonical incremental construction approach is needed to build 
only those simplices intersected by the splitting plane a. Instead of using a single list of active 
faces (AFL), the algorithm uses three disjoint lists containing: 

• AFL a : the (d — l)-faces intersected by plane a; 

• AFL i : the (d — l)-faces with all of the vertices in Pi ; 


AFL 2 : the (d — l)-faces with all of the vertices in P2', 







For each simplex a, the algorithm inserts its (d — l)-faces in the suitable face list. It then 
extracts faces (on which the next simplices will be built) from the AFL a alone; this ensures 
that each simplex built is part of the simplex wall 

The simplex wall construction process terminates when AFL a is empty. This process 
returns both T, a and the pair of active face lists AF L\ and AFL2- DeWall is then recursively 
applied to the pairs (Pi, AF L\) and {P2, AFL2), unless all the active face lists are empty. The 
splitting plane a is cyclically selected as a plane orthogonal to the axes of the E d space (X, Y 
or Z in P 3 ), in order to recursively partition the space with a regular pattern. Two-dimensional 
examples of the simplex wall construction and of the recursive application of the algorithm are 
shown in Figures 4 and 5, respectively. 

3.2 Uniform grid 

The DeWall algorithm is simple and easy to implement although in its naive implementation 
the asymptotic time complexity is not optimal nor is its practical efficiency good. An analysis 
of the algorithm shows that the main inefficiency is in the MakeSimplex function. 

Each simplex is constructed from an adjacent simplex face, by finding the dd-nearest point 
(i.e. the nearest according to the dd metric). This search entails performing an O(n) test for 
each simplex, where n is the number of sites in P. However, the construction of a new simplex 
in expected constant time is possible. 

The concept of local processing is often adopted in computer graphics either to speed up 
sequential algorithms or to achieve parallelism. The speed up technique proposed here is based 
on the E d extension of the uniform grid (UG) [13]; for simplicity, the use of the UG is described 
here for the case of DT in E 3 , supporting a regular partition of the space into hexahedral cells: 

UG={c ijk y, i,j,ke[0..N] (1) 

The main reason why uniform grid techniques are effective in geometric computations is that 
two points, which are far apart, generally have little or no effect on each other. A large class of 
geometric algorithms possess this property, ranging from visibility, to modeling (boolean oper¬ 
ations, intersection detection, etc.) and computational geometry (point location, triangulation, 
etc.) [14], 

The uniform grid is used as an indexing scheme for the fast detection of the dd-nearest 
point. A similar technique was also used by Fang and Piegl [15, 16] to speedup incremental 2D 








and 3D Delaunay triangulation. 

The space E 3 is partitioned into cubic cells following a regular pattern. The UG structure 
is built in a preprocessing phase, by computing for each cell c^jk the subset of points in P 
contained in c^k- 

The MakeSimplex function is designed such that, analogously to Maus’s proposal [17], the 
UG is scanned in order of increasing distance from /. Given this partial ordering of the sites, 
not all the points in P have to be analyzed for each face /. In fact, given a point p\ such 
that dd(f,pi) = di, all the points which are not contained in the sphere around / and pi 
will certainly have a dd value greater than d i, and it is pointless to evaluate their dd value. 
The analysis of the cells of UG can be stopped when there are no more cells contained in the 
circumsphere around / and the current dd-nearest point (Figure 6). 

The cells scanning order used is simpler than that proposed by Maus. Indeed we do not 
test the cells contained in circumspheres with increasing radius (the sphere_to_cells conversion 
is not a simple task) but we simply select and test all of the cells contained in the smallest 
cube circumscribed to each circumsphere. This method is simpler because it avoid the scan- 
convertion of spheres, and the number of cells selected is not much higher. Note that if the sphere 
radius selected is small (up to three times the cell edge length) the discretized circumsphere 
and the circum-cube are identical. 

The choice of the right resolution for the uniform grid space crucially affects the efficiency 
of the algorithm. In the reported implementation, the resolution of the UG is defined such that 
the number of cells is equal to the number of sites. 

3.3 DeWall Time Complexity 

The worst-case time complexity of the DeWall algorithm may be misleading: both the two 
techniques used (D&C strategy and Uniform Grid optimization) does not guarantee worst 
case optimality while offering good performances in pratical situations. It is possible to define 
patological datasets which cancel the efficiency of both the D&C strategy and the UG: if DeWall 
is applied to the dataset depicted in Figure 7, the construction of the first wall originates 
the entire triangulation (all the simplices in the triangulation intersect the splitting plane a); 
analogously, it is possible to choose site distributions that make the Uniform Grid not useful 
at all. In these pathological situations the DeWall algorithm reduces itself to an incremental 
construction algorithm, yielding a 0(n^~^ +1 ) worst case time complexity. In spite of this 
result, the algorithm behaves well in pratical cases (as shown in Section 4) yielding, in the 








three-dimensional case, a plain subquadratic behaviour versus a 0(n 3 ) worst case complexity. 

3.4 DeWall Space Complexity 

The algorithm space requirements are bounded by the space complexity of: 

• the point set P ; 

• the active face list AFL ; each AFL(n, d) is always a set of connected (d— l)-faces forming 
a unique (d — 1) surface in E d . Recalling that the number of (d — l)-faces of a polytope 
in E d of n vertices is at most 0(n^i J) [18], the worst case space complexity of AFL(n, d) 
is 0(nLlJ) ; 

• the outcoming triangulation; however, like the incremental construction algorithms, De- 
Wall can return each simplex as soon as it is built, avoiding explicitely storing the trian¬ 
gulation at run time. 

Therefore, the worst case space complexity of DeWall is 0(n^-~ J). The worst case size of 
the triangulation E in E d is 0(n^ “1) ; so it is interesting to note that the maximum space 
required by the algorithm in this worst case is lower than (or at most equal to in E 2 ) the 
size of the outcoming triangulation. On the other hand, on line triangulators need the current 
triangulation to be stored which is generally represented by the use of a hierarchical structure 
which holds the history of the construction process for fast point-in-triangle computations. 

4 Results and empirical evaluation 

The performance of the algorithm was tested on two classes of datasets. The first class consists 
of uniform datasets, where the locations of sites are generated using a uniform probability 
distribution function (Figure 8). In the second dataset class, the sites are organized into a 
number of bubbles with the density of sites decreasing as the distance from the bubble center 
increases (Figure 8). The sites in each bubble are generated using an approximation of a normal 
probability distribution function. 

For each dataset class and for each resolution (number of sites), a number of different 
datasets were generated in E 3 ; the times reported in Tables 1 and 2 are the means of the 
run times measured on each dataset. The machine used for the timings was an SGI Indigo 
workstation (MIPS R4000 cpu); the times include the uniform grid preprocessing. The results 
obtained show an empirically estimated complexity which is clearly subquadratic in E 3 . 
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Another way to empirically evaluate DeWall is to compare it with other implementations. 
We tested DeWall against two efficient Delaunay triangulators that are publicly available: 

• Incode: a totally incremental construction algorithm, with and without the use of the 
UG optimization technique 1 . Incode was implemented by using most of the DeWall’s 

• Qhull: a general dimension code for computing convex hulls and Delaunay triangulations. 
It is an implementation of the Quickhull algorithm [19] for computing the convex hull 2 . 
It was chosen because it qualifies as the fastest convex hull code for large datasets defined 
in low dimension spaces; 

• Detri: as part of the alpha-shape software, Detri builds the 3D DT by adopting an 
incremental insertion and flip approach [7] 3 . 

The results in Tables 1 and their graphical representation in Figure 9 show that DeWall 
is the most efficient of the four software programs on regularly distributed datasets, while it 
gives slightly slower times than Qhull on the bubble datasets. This is justified by the lower 
speedup obtained by adopting a UG on irregularly distributed datasets; the bubble datasets 
contains the worst distribution of sites for algorithms that use a UG (and therefore the DeWall 
algorithm). 

Some statistics on the execution of the DeWall algorithm on the uniform dataset are also 
reported in Table 1. The total number of tetrahedra returned is considerably lower than the 
theoretical upper bound in E 3 , 0(n 2 ): it was linear with the number of sites (approximately 
7 * n) in our experiments. The growth of the number of tetrahedra in the first wall is clearly 
sublinear (approximately 0(n~)). 

The mean number of cells visited for the construction of each simplex is not constant but 
shows a low increase with the dataset resolution. This is due to the fact that, for each face / on 
the Convex Hull (P) all of the cells contained in the positive halfspace of / have to be tested. 

1 Incode and DeWall are available in public domain at the address 

http://miles.cnuce.cnr.it/cg/swOnTheWeb.html 

Qhull is provided by the Geometry Center, University of Minnesota; the Qhull software may be 
downloaded from the WWW site http://freeabel.geom.umn.edu/software/download/qhull.html 

3 Detri is provided by the Software Development Group at the National Center for Supercomputing 
Applications (NCSA); info may be downloaded from the WWW site 

http://www.ncsa.uiuc.edu/SDG/Software/Brochure/Overview/ALVIS.overview.html 
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Figure 1: Merging of two partial DT in E~ space. 
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Figure 2: An example of DT in E 2 : a is the dividing line, and (the set of gray triangles) 
is the associated simplex wall; Si and S 2 are the triangulations returned by the recursive 
invocation of the DeWall algorithm on the two point set partitions. 



Function DeWall (P : point_set, AFL : (d-l)faceJList) : d-simplexJ_ist; 
var f : (d-l)face; AFL a , AFLi, AFL 2 : (d-l)faceJList; 

t : d-simplex; E : d-simplexJList; a : splitting_plane; 
begin 

AFL a , AFLi, AFL 2 : =emptylist; 

Point set _Part it ion (P, a, Pi, P 2 ); 

/* Simplex Wall Construction */ 
if AFL = 0 then 

t:=MakeFirstSimplex(P, a); 

AFL : = (d-1) faces (t) ; Insert (t,E); 
for each / G AFL do 

if Islntersected(f ,a) then Insert(f, AFL a ); 
if Vertices(f) C Pithen Insert (f , AFLi); 
if Vertices( f) C P 2 then Insert (f , AFL 2 ) ; 
while AFLa / 0 do 

f:=Extract(AFL a ); 
t:=HakeSimplex(f, P) ; 
if t / null then 

E:=E U {t}; 

for each /': /' G (d - 1 )faces{t) AND /' / / do 

if Islntersected(/',a) then Update (/',AFL a ) 
if Vertices( f') C Pi then Update (/', AFLi) 
if Vertices( f') C P 2 then Update (/',AFL 2 ) ; 

/* Recursive Triangulation */ 

if AFLi / 0 then E: =E U DeWall (Pi, AFLi) ; 

if AFL 2 / 0 then E: =E U DeWall(P 2 ,AFL 2 ) ; 

DeWall:=E; 

Procedure Update (f :face, L : face_list) ; 
begin; 

if Member(f,L) then Delete(f, L) 
else Insert (f , L) ; 

end; 


Figure 3: DeWall algorithm. 
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Figure 5: Some steps of the DeWall algorithm on a point set in E~. 
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